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1. Introduction 

Many gravity wave studies have been motivated entirely or in part to 
determine the perturbation introduced by the wave into the Richardson number 
field (e.g., Keller et al., 1983; Wurtele and Sharman, 1983; Thorpe, 1981; 
Fritts, 1979; Gossard and Hooke, 1975; Scorer, 1969). However, there have 
been few attempts at a systematic investigation of this problem, and since the 
Richardson number is a rather complex combination of the basic variables, it 
is difficult, even in the linear case, to form a conception of the Richardson 
number field, given the fields of the basic variables. 

Gossard and Hooke (1975, Section 36) derive the relevant formulas, and 
present a single calculation of a three layer model; but this model has 
infinite initial Richardson number and consequently is not well suited to 
compute Richardson number changes induced by gravity waves. 

In non-Boussinesq models, as is well known (Lamb, 1932), perturbation 
quantities increase in magnitude exponentially with height; and this effect 
has occasioned interest as a mechanism in the middle atmosphere for producing 
negative static stability (Hodges, 1967; Lindzen, 1981). Here we are 
concerned with the Boussinesq dynamics of the resonance lee wave, which dies 
out exponentially above its level of trapping. Since the trapping mechanism 
is wind shear, the upstream Richardson number is necessarily finite; however, 
the plane wave solution is not applicable. 

We propose here (1) to derive analytic solutions for an appropriate linear 
model with finite initial Richardson number, so that we may study the perturbed 
Richardson number field as a function of the mean flow parameters; (2) to 
execute numerical time-dependent simulations to compare with the analytic 
results; and (3) to obtain some useful qualitative non-linear formulations. 
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2. Linear analysis 

If ^ 1 is the perturbation stream function, w'the vertical velocity, u Q . and u' 
the mean and perturbation horizontal velocities, respectively, 6' the vertical 
displacement, and c' the vorticity, then we have, under steady state 
conditions, 


C* = vV 

w> “lx"* u o (z) I T" 

where subscript zpro indicates mean-flow quantities. 

The Boussinesq bouyancy, a' = 9 p'/p 00 > where p'is the perturbation 
density, and p Q0 a reference density, is related to these quantities by 


( 1 ) 


= N n 6 ' 
0 


■ V > 0 


2 

where N Q is the Boussinesq Brunt frequency, N 0 = -gdp 0 /p 0Q dz , for the mean 
density. 

We denote the perturbation stream function wave amplitude by 

A 

ip* ■= ip (z ) sin k x 


where k is the horizontal wavenumber of the disturbance, and similarly for other 
variables. The familiar vertical structure equation for the stream function 
is (Scorer, 1949): 




+ ( l 2 - k 2 )i = 0 , 
dz^ 


(3a) 
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where 




2 . N o 1 d2 % (3b) 

' U o^ 

The Richardson number for the perturbed flow may be written as 

R . yy/~ (tS/di i _ n_ I - dr’/tffc 

( 'iHf -i- ial )^~ f / +■ (chi! Jii)(d\i a lii) J ^ 

Mi di / 

2 2 

where Ri Q - N 0 /(du 0 /dz) is the Richardson number of the mean flow. In terms 

A 

of the vertical displacement 6 , (4) becomes, using (1), (2) and (3): 


n /n . _ i - t (d a /J b) + mVA/J-d?) it 

r\Jy 0 " “ A 2 

fit- J 2 m . 0 (JUcJdi-) 


(4a) 


A similar expression in terms of \p is 


.2, >i. 


Rl/Ki^ - — ITT . (4b) 

[ ) f i (iujAi) 


Equation (4b) reduces to equation (36-12) of Gossard and Hooke (1975) in the 
2 

case N q = constant. 

We now require a simple solution for a stable sheared flow permitting 
trapped waves, in which we may exhibit the pattern of the Richardson number 
field along side of the fields of \i>,‘ w, ui anda'. One candidate model is that 
of Palm and Foldvik (1960) for a flow in which the mean velocity is 
exponentially increasing with height. In this model, however, the mean 
Richardson number is therefore also exponentially varying, with the result 
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that it is difficult to interpret the perturbed Richardson number field. 
Interpretation is greatly facilitated if the basic Ri-field is uniform with 
height. 

An analytic solution for such a model has been presented by Wurtele 
(1953). This model assumes constant mean wind shear and stability, and 
therefore a constant Richardson number; however, it includes certain 
non-Boussinesq terms, for reasons to become apparent below. It has been known 
since the beginning of lee-wave theory (e.g., Scorer, 1949) that if both 
density gradients and potential termperature gradients are taken into account, 
and if 


p o (2) ■ p oo exp( ‘ Bz) 
and 

e o (z) = 9 oo exp(N o z/g) 
then the density-weighted variable 


satisfies the equation 


6 = const. 

2 

N = const, 
o 



( 6 ) 


Some now-familar approximations are involed in deriving this equation (Scorer, 
1949). The introduction of a constant mean wind shear into the model, 

u Q = U ( 1 + z/L ) U, L constant 
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thus provides a constant mean Richardson number 

Ri = N 2 L 2 / U 2 
o o 

against which one may study the perturbations arising out of the gravity-wave 
disturbance. If variables are rendered dimensionless by the scalings 

Z = 1 + z/L, <J>* = <p/U L , k* = k L , 6* = 6 L 

the vertical structure equation (6) becomes 


? * 

+ 4<t>* =0 
dZ^ 


( 7 ) 


where 


.2 , Rl 'o ..2 


4 ( 2 ) - 


- K 


and 


< 2 = k 2 + e|/ 4 , 


and equations (4a, b) for the Richardson number become 


!U,/fU c - 


1 _ ix 

[ 1 + '(l l kxj 


(8a) 


where 


and 


<5. = 


Ra/ R si a — 


?o \ l 6 


foe 


i 


Kj fj 3 ^ 


(8b) 
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The most convenient representation of the Richardson number field is the 
percentual change from the initial or mean field Ri, 


ARi _Ri_ 

R?“ ~ Ri n 
o o 


1 . 


From equations (7), we see that this quantity must be in phase with the fields 

of ‘ and u,' and 90 degrees out of phase with w. However, Ri is not linearly 

related to these fields, owing to the modification by the denominator. To 

see the effect of this shear term, consider the levels at which displacements 

d 6* 

have a maximum or minimum, i.e., — = 0. Then expanding the denominator by 

d Z 

the approximation 

(3+e)' 2 =l-2e + £ 2 , 


we have from (8a) 

= -2 i 2 6* sin kx + 3( i 2 6*sin kx ) 2 . 

Kl o 

Thus the percentual change in the Richardson number will be larger in 
magnitude when it is positive than when it is negative. We shall see in the 
simulations that this effect is important even for very small perturbations. 

i 

H £ 

The term in the numerator of (8a) is the linearized representation 

of the vertical stretching or shrinking (and hence the modification of the 

static stability) of an infinitesimal column of air when vertically displaced. 

2 

The term £ <5* in the denominator obviously arises from the structural equation 

(7), which gives the perturbation shear in terms of stream function or vertical 

velocity. The factor Ri 0 in this term suggest that it will dominate over the 
2 

modification of N in the numerator for the normal range of atmospheric 
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Richardson numbers. It also follows that the effect of the shear decreases 
with elevation and wave number. As will be shown, this hypothesis is borne 
out quantitatively in both analytic solutions and simulations. 

The solution to equation (6) regular as Z 00 may be written simply as 
(Wurtele, 1953) 


*. * +00 ^ V* 2 ’ 1 (9) 

o 

where y = Ri Q - 1/4, and is a modified Bessel function of imaginary 
order. This function is oscillatory in its argument for small argument and 
falls off exponentially for large argument. It is discussed and diagrammed in 
Appendix A. 

Thus, depending on the relative magnitudes of v» and k, the solution (9) 
permits a finite number, possibly zero, of gravity wave free modes, forced at 
z-r 0 ( Z = 1) and trapped at elevations inversely proportional to their wave 

lengths, and dying out exponentially above their level of trapping. The free 
modes (at z = 0) are given by 


K. ( 


) = 0 


( 10 ) 


that is, solutions satisfying a homogenous condition at the boundary. We 
shall be concerned only with conditions for which one or two wave solutions, 
or free modes, exist. These are determined by evaluation of the Bessel 
function by computer, as outlined in Appendix A. 
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The effect of the non-Boussinesq term 6 now becomes evident. The 
function K^(k) is singular as<-»- 0, becoming infinitely oscillatory (see 
Appendix A). In this model, regardless of how small k is, k will be limited 
below by 0/2. Thus the non-Boussinesq effect is to filter out all waves 
longer than the gravity wave asssociated with the mean density gradient. 
Another instance in which this effect permits a mathematically consistent 
solution occurs in the internal wave study by Mowbray and Rarity (1967). Our 
simulation model is entirely Boussinesq in its formulation; but in a 
simulation, of course, the longer waves (k 0) are eliminated simply by 
virtue of the limited horizontal and vertical extent of the computational 
grid. Since the wave length associated with the wave number 0/2 is of the 
order of 100 km and greater, this constitutes no limitation to the present 
study. 

The wave solution is related to its kinematic forcing at the boundary by 
a technique that is too familiar to warrant repetition here (Wurtele, 1953; 
Palm and Foldvik, 1960; Queney, 1960). We assume a surface of the form known 
as the "Witch of Agnesi" 

h(x) = ^-2 
a +x 

A 

which has the convenient Fourier transform htk) = Ha exp ( -a | k | ). 

The integral expression for the perturbation stream function is then 

kiu QjV d t in) 

kj (*) 

We follow the procedure referred to above, using the method of residues, and 
taking the principal value of the integral when poles exist in such a way as 
to cancel the upstream wave. 


f '(*■!*) = ^ ^ 


kx- 


~oo 


- 8 - 



The result is a wave solution, provided (10) is satisfied for some k = k^ , that 
has the form 


i|j'(x,z) = -2 it a e sin kx 

L K i' u <‘l> 


(12) 


where K. is the derivative of K. with respect to its argument. 
ip 


If other roots of (10) exist, these are of the same form and the 
solutions are additive. There is of course, in addition, a forced monotonic 
part of the solution that dies out exponentially both upstream and downstream, 
and does not concern us here. 


In the linear analytic computations we have calculated i|>' from (12), 
other variables from the linear relationships, and the Richardson number 
field from (8). 

The comparison numerical simulations are computed not for purposes of 
verification, but in order to identify limits on the assumption of linearity. 
As will appear below, this is not entirely a straightforward matter. The 
simulations are computed with the Boussinesq code described by Sharman and 
Wurtele (1983) and more fully by Sharman (1981). In every case, the 
computations are performed with a linear surface boundary condition, but with 
all non-linear terms present in the dynamic equations. The comparisons with 
the analytic computations thus isolate the departures from linearity due to 
the dynamics from those inherent in the calculation of the Richardson number 
itself. 
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3. Results of the analysis and simulations. 


Linearity in lee-wave problems is conventionally assessed by the parameter 
NH/U and we accept this criterion for the moment. For a sensitivity study 
of this parameter under the conditions N = const, U = const, see Miles (1965). 
Results will be presented for a variety of cases which are defined by Ri Q and 
NH/U. (See Table I for summary of cases and figures.) 

Figure 1 (a,b,c,d) presents the analytically derived perturbation fields 
for the conditions of Case I described in the legend. Note that as mentioned 
above in this and all computations herein, the Brunt frequency and wind shear 
are uniform, so that the basic Richardson number is also uniform with height. 

The conditions assumed for Case I specify Ri 0 = 8 and NH/U = 0.1, a 
highly linear problem, as will be seen below. Figures la and lb show, 
respectively, the perturbation stream function and the vertical velocity. 
Figure lc contains contours of total horizontal velocity, and the smallness of 
the perturbation relative to the mean flow is evident. In the areas near the 
surface where the total horizontal speed is reduced (i.e., where the contours 
of Figure lc are lifted), the shear is increased, and the Richardson number is 
decreased. Conversely, where the total horizontal speed is increased, the 
shear is decreased, and the Richardson number is increased. This qualitative 
assessment is borne out in Figure Id, which presents contours of the 
Richardson number. Quantitatively, however, we see, consistent with the 
analysis of Section 2, that even in this highly linear example, the fields of 
increased Ri and decreased Ri are not symmetric, the maximum increase being 
about 4.1 and the maximum decrease about -2.6. The trapped gravity lee wave 
is more stabilizing than destablizing to the basic flow. 
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For comparison, the simulation with nonlinear terms acting is presented in 
Figure 2, for Case I. By 600 time steps (9000 sec) a steady state has been 
achieved in the grid, the computational radiation boundary condition (Sharman 
and Wurtele, 1983) permitting the wave to be established all the way to the 
outflow boundary without significant reflection. Figure 2a shows the 
streamfunction, here with the phase determined by the position of the 
lOOm-obstacle. The wavelength (about 20 km) in the simulation agrees with the 
analytic solution of Figure la. Figure 2a shows a similar agreement of the 
vertical velocity fields, where the analytic solution has a maximum of 0.4 m/s 
and the simulation, 0.5 m/s, both attained at an elevation of 3 km. Whether 
this is scaled by the surface mean wind speed of 10 m/s or the mean speed at 
3 km of about 21 m/s, the perturbation must be considered small. However, in 
Figure 2c, corresponding precisely to Figure Id, we see the same lack of 
symmetry in the Ri-field exhibited by the analytic solution. The simulation 
Richardson number changes are about the same as those obtained by analysis. 
Here the greatest increase in Ri is about 2.8 and the greatest decrease is 
about -1.6. 

We may now consider Case II, for which N and U are the same as in Case I, 
but H = 500 m. Hence Ri Q = 8, as before, but NH/U = 0.5. Here one would 
expect to see the linear result of Case I amplified, by a factor of about 
five, but with no significant nonlinear effects. And this is, in fact, what 
we find, in all fields except that of the Richardson number. Figures 3a,b,c,d 
present, respectively, the stream function, vertical velocity, horizontal 
velocity and Ri-field for the analytic solution of Case II. The larger 
disturbance shows more clearly the close correspondance between the maximum 
shear (Fig. 3c) and the minimum Richardson number, and vice versa. For Case 
II, the analytic Ri-field is highly asymmetric and varies between 3.2 and 32. 
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Figures 4a,b,c show the simulation fields for Case II. Again the stream 
function field (Figure 4a) is in agreement with the linear analytic solution. 
However, the most convincing quantitative evidence of linearity is in the 
vertical motion fields of Figs lb, 2b, 3b and 4b. The forcing, as we have 
seen, increases from Case I to Case II by a factor of five. The linear 
solutions (Figs, lb, 3b), not surprisingly, show an increase of the maximum 
amplitude of vertical velocity from 0.4 to 2.0 m/s. The simulations 
(Figs. 2b, 4b) show a corresponding increase from 0.5 to 2.7 m/s. In 
contrast, the Richardson number field (Figure 4c) shows strong nonlinearity. 

The positive cells are much larger in size than the negative cells, and the 
field varies from a maximum of 111 to a minimum of about 1.5. The nonlinear 
dynamics of the simulation, although not evident in the other fields, has 
enhanced the asymmetry of the analytic solution for the Richardson number 
field of Fig. 3d. As we have noted, this enhancement did not occur in the 
highly linear Case I. 

Although the Richardson number of the basic flow determines the 

dimensionless wave length, other quantities, both dimensional and 

nondimensional (including dimensional wave length) will vary with stability and 

shear, even though the basic Richardson number is fixed. As an illustration of 

-2 -1 

this we consider Case III, for which N 0 = 0.5 x 10 s , L = 5.66 km and 

H = 500 m. Thus Ri = 8, as in Cases I and II, but both shear and stability have 

o 

been decreased by a factor of two. 

Figure 5 (a,b,c,d) displays the fields from the linear analysis for Case 
III. The dimensional wave length has increasd to about 42 km. The level of 
trapping has also increased; this is most evident in the vertical velocity 
field Figure 5b, where the maximum is centered at about 7 km elevation, instead 
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of at 3 km for Cases I and II. Since the level of trapping is higher, the 
energy from the forcing is supplied to a larger volume of atmosphere, and the 
disturbance amplitude is correspondingly smaller. 

So far we have discussed conditions permitting a single trapped wave. In 

order to limit the system to a single free mode, it was necessary to specify a 

relatively small Richardson number by atmospheric standards, Rio = 8. 

Reference to Appendix A and to Figure A1 show that doubling this value to 

Ri = 16 defines a condition in which two free modes are admissible, of 
o 

wavelengths approximately 15 km and 30 km. To define Case IV then, we take 
N = 10" 2 s _1 , L = 4 km, and H = 500 m. Thus Ri 0 = 16 and NH/U = 0.5. As with 
Case II, we shall see that all fields except that of the Richardson number 
exhibit linear behavior. 

Figure 6 (a,b,c,d) displays the fields corresponding to those of Figures 
1,3 and 5, that is, the stream function, vertical velocity, total horizontal 
velocity, and total Richardson number for the first mode only. The patterns 
are the familiar ones from The Ri Q = 8 cases. And the same is true of Figure 6 
(a,b,c,d) containing the corresponding representations for the second mode 
only. The second mode, of course, has two cells in the vertical in the 
streamfunction, vertical velocity, and Richardson number fields, and three 
cells in the horizontal velocity field. This mode has its maximum amplitude 
at a higher level than the first mode velocity fields, but the maximum are 
less than those of the first mode. Richardson number in each mode is reduced 
by slightly more than 50% in the centers of the negative cells. 

The total linear free-wave solution is the sum of these two modes, 
presented in Figure 8 (a,b,c). The perturbation stream function field of 
Figure 8a shows clearly the superposition of a short and long wave, the former 
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dominating at lower, and the latter at higher elevations. The same is true of 
the vertical velocity field of Fig. 8b. Yet the wavelengths of the two waves - 
derived, it will be recalled, as roots of a Bessel function - are not of 
integral ratio, and thus the combined pattern is not strictly periodic. In 
the figure, there is a maximum value for w of about 2.3 m/s. Thus the 
problem is within the bounds of linearity, and the fields of these two 
figures have an ordered pattern. 

The Richardson number field (Fig. 8c), however, is complex and the pattern 
would be difficult to anticipate. The maximum Richardson number in the 
graphed portion of the field is of the order of 10®, but the minimum, 3.1, is 
greater than the minima of the two modes separately. 

4. Limits of linearity 

For the trapped lee wave, as the wind shear increases, the level of 

trapping is displaced downward, and wave amplitudes are increased. As the 

amplitude becomes larger, nonlinear effects may be expected to become more 

significant. As an illustration of this we consider Case V, for which we 
-2 -1 

choose N = 2 x 10 s ,L= 1.414 km and H ■ 500 m. Thus Ri Q = 8, as in Cases 
I and II, but both shear and stability have been increased by a factor of two. 
And NH/U = 1, also twice the value of Case II. The Figure 9 (a,b,c,d) shows 
clearly the increased magnitude of all quantities. The maximum vertical 
velocity is now 32% of U and the minimum Richardson number is 0.51, 
approaching the critical. However, traditionally a value of unity for NH/U 
would be said to render a linear analysis invalid, and we must at this point 
consider the problem of criteria for nonlinearity. 


- 14 - 



We first note that the conventional parameter of linearity, NH/U, 
developed for uniform flow when these three parameters determine the flow 
completely, is less easy to justify when shear is present. Another length, L, 
is available, and H/L would seem to be a candidate dimensionless parameter of 
the problem. The most thorough perturbation analysis of shear flow is that 
developed by Brown and Stewartson (1982) for an essentially nonlinear problem, 
that of a gravity wave propagating toward a critical level. Although the 
forcing in their problem is not specified, if it is identified with kinematic 
forcing at a lower (or upper) boundary, their expansion parameter can be 
identified as our H/L. From equation (12) we have an expression for the 
vertical velocity w : 


w 

U 


= -2 ir (akj) e 


-aki 


H K 1u< *1 Z ) 
L 


cos 


kjX 


It is evident here that H/L is the ratio determining the amplitude, and that 
the stability, as represented by N q , plays little role. It is true that the 
quantity 

(akj) exp (-akj) , 


associated with the shape of the forcing obstacle, is a function of N q , but it 
is not a monotonic function, and its maximum value is e~* 


Thus in general we have the results that NH/U is the amplitude parameter 

u d u 0 

for the zero-shear case and H/L, or more generally, — -g— v is the 

o 

amplitude parameter for a shear flow. However, we cannot pretent that the 
problem has been solved. In the first place, the shear flow must depend on 
some combination of these parameters; otherwise L -► °° (shear approaching 
zero) gives a wrong result. The failure of the criterion for shear flow to 
pass into the corresponding result for zero-shear flow is characterstic of all 
theoretical formulations. 
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In the second place, what does "nonlinearity" mean in this context. In 
lee wave theory it has been taken to mean ii' -*■ -u Q , 9p/3z-*-0 , i.e., 
overturning and wavebreaking. In order to examine the consequences of such a 
situation, we have used a simulation code better adapted to extreme nonlinear 
parameter- space, that of Pihos and Wurtele (1981). In Figures 10-12a, b, we 
present the streamline pattern and density contours for Ri Q = 8, under 
conditions of Case V, except that NH/U takes on the values of 1.0, 1.5, and 
3.0, respectively. We see that (1) overturning - as evidenced by heavier 
fluid above lighter - is not present until NH/U = 3; and (2) even at this 
value, the instability is highly local and does not spread with increasing 
time. In fact, the choice of values to be plotted in the streamline field 
fails to resolve the overturning at all. 

For all these three runs, the value of L is 4 km, so that H/L is 0.174, 
0.265 and 0.530 respectively. The beginnings of instability for this last 
value suggests that H/L is the decisive dimensionless parameter for stratified 
shear flows. 

As a control, we submit the experiment of Figure 13, in which there is 
zero shear, i.e. U = constant, N = constant. In this case, as indicated 
above, the parameter NH/U governs linearity, and the value of unity assigned to 
it in this run has clearly produced overturning. 
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5. Nonlinear Analysis 

It is not clear that a linear (or perturbation) analysis is the best or 
only possible approach to this problem. Some light is cast by following a 
line begun by Scorer (1969), using one Lagrangian coordinate, say, potential 
temperature. The transformation equations - somewhat more complicated than in 
the hydrostatic model - are listed for convenience in Appendix B. We consider 
only the steady-state, in which isentropes are streamlines. We thus have from 
(B5), the equation of continuity, and (B6), the vorticity equation, in terms 
of the total (finite) variables: 


3 / 3z \ _ 

37^ u W ” 0 ’ 


(13) 


where z is the height of an isotropic surface or streamline, and 


..H = ifLLr 1 

u a x " e k 3 e ; 3x 


(14) 


where c is the vorticity. The first of these may be integrated to yield, 
using (B3), 


u_ 

u. 


(15) 


where U Q and N Q are values far upstream on the streamline, both functions of 
height. This equation simply quantifies our intuitive knowledge that when a 
streamline becomes vertical (u = 0), the stability also vanishes. This is 
not in general the case in linear analysis, if we identify u Q + u' with u 
and N 2 -3 o'/3z with N 2 . From (1) to (3) we obtain in contrast to (15), 
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u i|>* 3 / 

u m 2 3 z 1 

o N 



This equation illustrates the familiar result that 1 inear sol utions satisfy 

IT 

the nonlinear equations provided - — is constant. For the model of Section 

o 

2 , 





This will, in general, occur in a finite layer, with and N 2 both near 
zero or negative, emphasizing that when convective instability enters the 
problem, the Richardson number ceases to be a useful concept. 

We may gain a little further insight from the vorticity equation (B6), 
again using (B3) 


X = 
3x 


- N 2 
n 3x 


2 

Since by (15) the ratio N /u is constant along a streamline, this may be 
integrated to yield 


£u-c 0 > - -« 


where 6 is the vertical displacement of a streamline. Qualitatively this 
equation states that (in this right-handed system) the vorticity increases 
when the streamline is depressed and decreases when it is lifted; but it 
further states that the change is linear with streamline displacement and 
specifies the constant of proportionality for each streamline. 
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6. Conclusions 

An analytic solution has been presented for a stratified fluid of 
arbitrary constant Richardson number. By computer aided analysis the 
perturbation fields, including that of the Richardson number can be 
calculated. The results of the linear analytic model were compared with 
non-linear simulations, leading to the following conclusions. 

1. The perturbations in the Richardson number field, when small, are produced 
primarily by the perturbations of the shear. 

2. Perturbations of in the Richardson number field, even when small, are not 
symmetric, the increase being significantly larger than the decreases. The 
linear analytic solution and the nonlinear simulations both confirm this 
result. 

3. As the perturbations grow, this asymmetry increases, but more so in the 
nonlinear simulations than in the linear analysis. 

4. For large perturbations of the shear flow, the static stability, as 

2 

represented by N , is the dominating mechanism, becoming zero or negative, 
and producing convective overturning. 

5. The convectional measure of linearity in lee wave theory, NH/U, is no 

H du 0 

longer the critical parameter. It is suggested that g — gp takes on, this 

o 

role in a shearing flow. 
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TABLE I 


Case 

Figure 

Hs) 

U(m/s) 

I 

1,2 

.01 

10 

II 

3,4 

.01 

10 

III 

5 

.005 

10 

IV 

6,7,8 

.01 

10 

V 

9,10 

.02 

10 

VI 

11 

.02 

10 

VII 

12 

.02 

10 

VIII 

13 

.02 

10 


L(km) 

R l°_ 

H(km) 

NH/U 

H/L 

2.83 

8 

0.1 

.1 

.035 

2.83 

8 

0.5 

.5 

.177 

5.66 

8 

0.5 

.25 

.088 

4.00 

16 

0.5 

.5 

.125 

1.41 

8 

0.5 

1.0 

.354 

1.41 

8 

0.75 

1.50 

.265 

1.41 

8 

1.50 

3.00 

.530 

00 

00 

0.50 

1.00 

0 
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LEGENDS 


Figure la. 


Figure lb. 


Figure lc. 


Figure Id. 


Figure 2a. 
Figure 2b. 


Figure 2c. 


Figure 3a. 


Figure 3b. 


Figure 3c. 


Graph of perturbation stream function from linear analysis in 
Case I (see Table I). All units in this and other computations 
are MKS. 

Vertical velocity from analysis for Case I. Contour interval is 
0.2 m/s. In this and all subsequent figures of vertical velocity 
from analysis the first (leftmost) cell is negative. 

Total horizontal velocity from analysis for conditions of Case I. 
Contour interval is 2 m/s. 

Richardson number field from analysis for Case I contoured for 
the quantity (Ri - Ri 0 )/Ri Q . The contour interval (dimensionless) 
is 0.05. In this and all subsequent figures of Ri from analysis 
the first (leftmost) cell is a region of Ri increase. 

Stream function in simulation for conditions of Case I, at time 
step 600 (9000 sec). 

Vertical velocity in simulation of Case I. Contour interval is 
0.1 m/s. (Note that contour interval is one-half that of Fig. lb). 
Richardson number field in simulation of Fig. 2a. As in Fig. lc, 
the quantity contoured is (Ri -Ri 0 )/Ri Q , and the contour interval is 
0.05. 

Graph of perturbation stream function from linear analysis for 
Case II. (See Table I) . 

Vertical velocity field from linear analysis for Case II. Contour 
interval is 0.2 m/s . 

Total horizontal velocity contours from linear analysis for 
Case II. Contour interval is 2 m/s . 
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Figure 3d. 
Figure 4a. 
Figure 4b. 
Figure 4c. 

Figure 5a. 

Figure 5b. 

Figure 5c. 

Figure 5d. 

Figure 6a. 

Figure 6b. 

Figure 6c. 

Figure 6d. 

Figure 7a. 
Figure 7b. 


Total Richardson number field from linear analysis for Case II. 
Contour interval is 1.0. 

Graph of perturbation stream function from simulation for Case II, 
600 times steps (7200 sec). 

Vertical velocity contours from simulation for Case II. Contour 
interval is 0.2 m/s. 

Richardson number field from simulation for Case II. The quantity 
contoured is (Ri-Ri Q )/Ri o with a contour interval of 0.10, or 
equivalenctly total Ri with a contour interval of 0.8. 

Graph of perturbation stream function from linear analysis for 
Case III (see Table I). 

Vertical velocity contours from linear analysis for Case III. 
Contour interval is 0.2 m/s. 

Total horizontal velocity contours from linear analysis for 
Case III. Contour interval is 2.0 m/s. 

Total Richardson number field from linear analysis for Case III. 
Contour interval is 1.0. 

Graph of perturbation stream function from linear analysis for 
Case IV (see Table I). First mode only. 

Vertical velocity contours from linear analysis for Case IV, first 
mode only. Contour interval is 0.2 m/s. 

Total horizontal velocity contours from linear analysis for Case IV 
first mode only. The contour interval is 2 m/s. 

Total Richardson number field from linear analysis for Case IV, 
first mode only. The contour interval is 2.0. 

Same as Fig. 6a, but for second mode only. 

Same as Fig. 6b, but for second mode only. 
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Figure 7c. 
Figure 7d. 
Figure 8a. 

Figure 8b. 

Figure 8c. 

Figure 9a. 

Figure 9b. 

Figure 9c. 

Figure 9d. 

Figure 10a. 

Figure 10b. 
Figure 11a. 

Figure lib. 
Figure 12a. 


Same as Fig. 6c, but for second mode only. 

Same as Fig. 6d, but for second mode only. 

Graph of complete perturbate in stream function from linear 
analysis for Case IV. Sum of the two modes. 

Vertical velocity contours from linear analysis of two mode sum 
for Case IV. Contour interval is 0.2 m/s. 

Total Richardson number field from linear analysis of two mode sum 
for Case IV. The contour interval is 2.0. 

Graph of perturbation stream function from linear analysis for 
Case V (see Table I.) 

Vertical velocity contours from linear analysis for Case V. 

Contour interval is 0.2 m/s. 

Total horizontal velocity contours from linear analysis for Case V 
Contour interval is 2.0 m/s. 

Total Richardson number field from linear analysis for Case V. 
Contour interval is 1.0. 

Stream function at 800 time steps (4800 sec) as simulated by the 
non-linear gravity wave code of Pihos and Wurtele (1981) for 
Case V. 

Total density field corresponding to Figure 10a. 

Stream function at 600 time steps (3600 sec) as simulated by the 
non-linear gravity wave code of Pihos and Wurtele (1981) for 
Case IV (see Table I). 

Total density corresponding to Figure 11a. 

Stream function at 900 time steps (2700 sec) as simulated by the 
non-linear gravity wave code by Pihos and Wurtele (1981) for 
Case VII (see Table I). 
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Figure 

Figure 


12b. Total density field corresponding to Figure 12a. 

13. Stream function at 400 time steps (4000 sec) as simulated by the 
non-linear gravity wave code by Pihos and Wurtele (1981) for 
Case VIII (see Table I). 
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APPENDIX A 


The modified Bessel function K V (S) is discussed and tabulated for real 
v and £ in may sources (e.g. Abramowitz and Stegun, 1964; Sections 9. 6-9. 8). 
However, the useful function (5), although arising naturally in fluid 
dynamical contexts (e.g. Booker and Bretherton, 1967; Mowbray and Rarity, 

1967) seems to have been tabulated only by Morgan (1947). Morgan’s work is 
relatively difficult of access; and since it was done at a time when computing 
capability was relatively slight, his tables provide less than adequate 
resolution. Therefore, we present here some relevant formulas and a diagram 
illustrating the behavior of the function. 


A qualitiative characterisation is easily obtained by means of the 
integral representation 


fyd) = 


00 



O 


or/*,£) dX 


(Al) 


derived from Abramowitz and Stegun (1964, 9.6.24) by the substitution v = iy. 
From this formulation, by straightforward application of the method of 
steepest descents, we derive the following two asymptotic limits, for y»£ 
and 'P 5 . For P » 5 , 
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These are given by Morgan (1947). Thus in the range y^£ , (£) changes 

from oscillatory to exponential behavior as £ increases. The oscillatory 
singularity in (A2) for £-*-0 is familiar as arising in the problem of the 
propagation of a gravity wave toward a critical level (Booker and Bretherton, 
1967). 

The values of K^(C) for use in Section 2 were calculated using three 

different relations for various ranges of the parameters y and £ . For large £ 

2 r 2 

(|£| > 20, or | £ | > 5 and 2 1 £ | > y , or |^-| >5y), the asymptotic relation (A3) 
with higher order terms included (see e.g. Abramowitz and Stegun Eq. 9.7.2) 
was used. For large y ( £>20, y>£ ) Debye’s asymptotic expansion (A2) with 
higher order terms included (see Abramowitz and Stegun Eq. 9.3.7) gave 
satisfactory results. For all other parameter values the standard ascending 
series (Abramowitz and Stegun Eq 9.1.10) was used to compute J. (i£) , with the 
value of K iy (£) provided by the identity 

If) ‘ -*5&J f «m) (■<)] 

The leading term of the result of this process reduces to 

fry. W = {phCf*) “*(/*** l- vyf'fy)] 

where r represents the gamma function. This representation is also given in 
Morgan (1947). 

The most convenient form for visual representation of the function K (£) 

•iy 

is a contour plot in (y,£)-space, following Jahnke, Emde, and Losch (1960, 

Figs. 79, 87) for Bessel functions of the first kind and Neumann functions. 

This plot is shown in Fig. Al. The figure readily shows the change from 
oscillatory solutions to exponentially decaying solutions in the neighborhood 
of y = £. In the oscillatory region, for the given order ( = /Ri Q in our 
model), the function becomes infinitely oscillatory as the argument £ 0. 
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LEGEND FOR APPENDIX A 


Figure Al. Contour plot of the Bessel function K^(£) in the ( E, , y ) 
plane. Negative regions are shaded. The format is selected for compatibility 
with Figures 79 and 87 of Jahnke, Emde, and Losch (1960). 
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APPENDIX B 

A transformation to (x,6,t) coordinates requires the following identities: 



All variables are total quantities and not perturbations. 

^ For completeness we list the momentum equations In (x,6)-coord1nates 
although they will not be used here: 



where it is the barotropic pressure function 

C P (PjPc) ^ 
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The continuity equation in the non-hydrostatic system becomes time-dependent. 


I /ii I j. ^ / sCC ii ) - 0 


(B5) 


if the vorticity is defined as 


C 


3u~ _ £<c 
3H )i 


the vorticity equation 


11 




v. \7 £ = 


becomes in the ( x,e)-system 


± b£ 
6 >/■ 



